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Summary 

This paper presents a method for calculating vis- 
cous effects on two- and three-dimensional unsteady 
transonic flow fields. An integral boundary-layer 
method for turbulent viscous flow is coupled with 
the transonic small-disturbance potential equation in 
a quasi-steady manner. The boundary- layer calcula- 
tion uses Green’s lag-entrainment equations for at- 
tached flow and an inverse boundary-layer method 
for flows with mild separation. Three-dimensional 
viscous effects are approximated by a stripwise appli- 
cation of the two-dimensional boundary-layer equa- 
tions. The method is demonstrated for several test 
cases, including two-dimensional airfoils and a three- 
dimensional wing configuration. The applications for 
two-dimensional airfoils include an example that il- 
lustrates the method for calculating aileron buzz and 
thus demonstrates the present method for analyzing 
a key aeroelastic problem. Comparisons with invis- 
cid calculations, other viscous calculation methods, 
and experimental data are presented. The results 
demonstrate that the present technique can econom- 
ically and accurately calculate unsteady transonic 
flow fields having viscous- inviscid interactions with 
mild flow separation. 

Introduction 

Computational methods for accurately calculat- 
ing unsteady transonic flow for aeroelastic applica- 
tions are rapidly maturing (ref. 1). For example, 
Malone, Sankar, and Sotomayer (ref. 2) calculated 
unsteady air loads on the F-5 fighter wing with a full- 
potential computer code, Steger and Bailey (ref. 3) 
calculated aileron buzz with a Navier-Stokes code, 
and Anderson and Batina (ref. 4) calculated un- 
steady pressure distributions for both two- and three- 
dimensional configurations with an Euler code and 
a transonic small-disturbance potential code called 
CAP-TSD (Computational Aeroclasticity Program- 
Transonic Small Disturbance) (ref. 5). Other appli- 
cations of Euler codes and Navier-Stokes codes il- 
lustrate the complex flow phenomena that can be 
computed by these methods. However, full-potential, 
Euler, and Navier-Stokes codes usually require large 
amounts of computer time and, as a result, are cur- 
rently too expensive for routine applications. Thus, 
substantial efforts have been devoted to the develop- 
ment of transonic small-disturbance codes (ref. 6). 

For flows involving weak or moderately strong 
embedded shock waves, inviscid calculations that 
use the TSD equations have produced accurate so- 
lutions in many cases: thin airfoils (ref. 7), thin 

wings (ref. 8), wung-canard combinations (ref. 9), and 
realistic aircraft configurations (ref. 6). As shock 


waves increase in strength and move aft on the air- 
foil, viscous effects become significant and must be 
accounted for in the computations to obtain ac- 
curate solutions (ref. 10). For flows that remain 
attached, integral boundary-layer methods may be 
coupled with the inviscid analysis by viscous-inviscid 
iteration. These interactive boundary-layer tech- 
niques have produced viscous solutions that agree 
well with experimental results (refs. 11 to 14). 

For separated flows, integral techniques are also 
available. In particular, LeBalleur (ref. 15) devel- 
oped a fully unsteady viscous-inviscid integral tech- 
nique. Good results were achieved w^hen LeBalleur 
and Girodroux-Lavigne (ref. 16) applied the tech- 
nique to several airfoils that had strong viscous- 
inviscid interactions and extensive regions of flow 
separation. (See ref. 16.) The technique, however, 
can require large computer resources; some cases in 
reference 16 required up to 15 viscous-inviscid itera- 
tions at each time step to obtain converged solutions. 
Melnik and Brook (ref. 17) specialized LeBalleur’s 
technique, with some modifications, to steady cal- 
culations for inclusion in the GRUMFOIL computer 
code. Calculations made with this code agree rea- 
sonably with experimental data up to and slightly 
beyond maximum lift. 

This paper presents an efficient method for calcu- 
lating viscous effects on two- and three-dimensional 
configurations for unsteady transonic turbulent 
flow r s. The inverse boundary-layer method in refer- 
ence 17 is incorporated into the CAP-TSD computer 
code (refs. 5, 18, and 19) in a quasi-steady manner. 
Carter’s method (ref. 20) is used to couple the inverse 
calculations with the inviscid algorithm. Green’s lag- 
entrainment equations are included to calculate at- 
tached flows. The resulting computer code is applied 
to several test cases, including both two-dimensional 
airfoils and a three-dimensional wing configuration. 
The results demonstrate that the present technique 
can economically and accurately calculate unsteady 
transonic flow fields involving viscous-inviscid inter- 
actions with mild flow separation. 


Symbols and Abbreviations 

CAP-TSD Computational Aeroelasticity 
Program-Transonic Small 
Disturbance 


C E 


entrainment coefficient 


Cj skin-friction coefficient 

Cp pressure coefficient 


Cp 

normalized unsteady pressure 
coefficient; first harmonic 
of C p divided by oscillation 
amplitude 

C T 

shear stress coefficient 

C 

airfoil chord, m 

c l 

lift coefficient 

c m 

pitching-moment coefficient 
about quarter- chord 

EXP 

experimental 

f 

oscillation frequency, Hz 

hi • • • i h 

functions in transonic small- 
disturbance equation defined 
by equation (2) 

H, H u H 

boundary-layer shape factors 

k 

reduced frequency, ujc/2U 

M 

free-stream Mach number 

m 

= PeU e 6* 

JVpr.t 

turbulent Prandtl number 

N Re 

Reynolds number, Ucjv 


Sutherland number 

Q 

dynamic pressure, psf 

R\, i?2, #3 

constants in equations (9) 
to (11) 

r 

= iV 1/3 
iV PM 

S 

airfoil surface function 

t 

nondimensional time, 

At 

nondimensional time step 

t 

time, sec 

U 

free-stream velocity, m/sec 

Urn 

magnitude of reverse flow in 
boundary layer 

u v 

streamwise velocity of viscous 
flow in boundary layer 

x, y, z 

Cartesian coordinates in 
streamwise, spanwise, and 
vertical directions 

a 

angle of attack, deg 

«0 

mean angle of attack, deg 

ai 

dynamic pitch angle, deg 


a 

= 6*/6 

P 

flap angle, deg 

1 

ratio of specific heats 

A(...) 

indicates jump in . . . 

6 

boundary-layer thickness, m 

6* 

boundary-layer displacement 
thickness, m 

V 

= zjb 

V 

= 0 - »7*)/(l - V*) 

Ps 

fraction of semispan 

V* 

height of reversed flow region 

0 

boundary- layer momentum 
thickness, m 

V 

kinematic viscosity, m 2 /sec 

P 

density 


inviscid-disturbance velocity 
potential 

uj 

relaxation factor 

Subscripts: 


B 

aileron buzz 

€ 

boundary-layer edge 

i 

inviscid quantity 

le 

leading edge 

te 

trailing edge 

V 

viscous quantity 


All angles are positive for trailing edge down. 
Moments are positive for leading edge up. Hinge 
moments are taken about the hinge axis. 

Governing Equations 

The inviscid flow code used in this analysis is 
the transonic small-disturbance potential computer 
code CAP-TSD developed at NASA Langley Re- 
search Center by Batina et al. (refs. 5, 18, and 19). 
The CAP-TSD code uses an approximate factoriza- 
tion algorithm (ref. 18) for time-accurate solution of 
the unsteady TSD equation. The code has been ap- 
plied extensively to airfoils (refs. 4 and 18), wings 
(ref. 21), wing-body configurations (ref. 5), and com- 
plete aircraft configurations (ref. 6). These refer- 
ences include comparisons with experiments as well 
as with other computer codes for computational fluid 
dynamics. 
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The viscous analysis presented in this paper in- 
teractively couples the CAP-TSD inviscid flow code 
with an integral boundary-layer technique to model 
turbulent viscous flow effects. The direct boundary- 
layer method for attached flow is based upon Green’s 
lag-entrainment equations and is a modified applica- 
tion of the method described in reference 14. The 
equations are repeated herein for completeness. The 
inverse boundary-layer equations are based upon the 
work of Melnik and Brook (ref. 17) and are in- 
cluded in the CAP-TSD computer code in a stripwise 
manner. 


Inviscid Equations 

The CAP-TSD computer code solves the modified 
transonic small-disturbance equation in conservative 
form 

dfo , dfi d /2 , d /3 


+ — = 0 

dt dx dy dz 


( 1 ) 


where 0 is the inviscid- disturbance velocity potential: 


fo — ~A(f>t — Bcj) x 

(2a) 

fl = E\4> x + Fl<f> 2 + G\<f)y 

(2b) 

/2 — 4>y + H\(f>x<i>y 

(2c) 

h = <f>z 

(2d) 


The coefficients A, B, and E\ are defined as 
A = M 2 B = 2 M 2 Ei = 1 - M 2 


Choices for the coefficients F \ , G\, and H\ depend 
upon the assumption used to derive the TSD equa- 
tions. In this paper, the two-dimensional calcula- 
tions are made with the following “NLR” coefficients 
(ref. 22): 


n = - ^ [ 3 - (2 - 7 )M 2 ] M 2 
Hi = -M 2 


For the three-dimensional calculations, the following 
“NASA Ames” coefficients (ref. 22) are used: 

= -\{l + l)M 2 

Gi = \( 7-3)M 2 
Hi = -( 7-1 )M 2 


Also, the CAP-TSD code incorporates modifications 
to the coefficients in equations (2); these modifi- 
cations were developed by Batina (ref. 19) to ap- 
proximate the effects of shock generated entropy or 
vorticity. 

The boundary conditions on the wing and wake 
are 


4>f = Sf + 

(xie <X< X t e, Z = 0±) 

( 3 ) 

O 

II 

<1 

(x > x te ; * = o*) 

( 4 ) 

A(0 X + 0t) = 0 

(x > Xte\ Z = 0*) 

( 5 ) 


where the superscript ± refers to the wing upper or 
lower surface, the function S(x,t) denotes the wing 
surface, and A(. . .) indicates a jump in the bracketed 
quantity across the wake. In the far field, nonreflect- 
ing boundary conditions similar to the ones devel- 
oped by Whitlow (ref. 23) are implemented in the 
CAP-TSD code. References 6 and 23 contain details 
of the derivation of those boundary conditions. 


Viscous Equations for Attached Flow 

The effect of a turbulent viscous boundary layer 
for attached flow is modeled in a quasi-steady man- 
ner by Green’s lag-entrainment equations as imple- 
mented in reference 14. References 14 and 24 present 
additional details. The boundary-layer equations for 
attached flow are 

~ = \ c > - ( h + 2 - m < ) <6 > 


rj dC E 

dx 


- ' {hThT ' " ''1 ‘ (eIf) 


EQ 


( 8 ) 


Equation (8) for the entrainment coefficient is taken 
from reference 12 and differs slightly from the equa- 
tion given in reference 24. The surface velocity gra- 
dient 4> X x is smoothed for numerical stability dur- 
ing the computations as discussed in reference 14. 
The subscript e in these equations refers to quan- 
tities at the boundary-layer edge, the subscript EQ 
denotes the equilibrium conditions, and the subscript 
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EQO denotes the equilibrium conditions in the absence of secondary influences on the turbulence struc- 
ture (ref. 24). The various parameters in these equations (be., Cf, F, H , Hu M e , CV, r, (^:^) E q’ and 


(CV)eqo) are defined in the appendix. 


Viscous Equations for Separated Flow 

In flow fields that contain regions of separation, the bound ary- layer equations are written in inverse form. 
Thus, these equations can be solved when given a specified streamwise variation of boundary-layer displacement 
thickness as represented by a perturbation mass flow parameter m = p e U e S*. The solution to the inverse 
equations (i.c., the viscous velocity at the edge of the boundary layer) is then used in a relaxation formula 
to update the displacement thickness and calculate a new value of m. This iterative process is repeated at 
each time step until convergence is achieved. This particular inverse form of the bound ary- layer equations was 
developed by Vatsa and Carter at the United Technologies Research Center, and it is completely compatible 
with Green’s original lag-entrainment equations in regions of attached flow. The inverse equations are 


- 1 dm | 

1 dU e m ~dx 

lid d!ii 1 

{Cz-lCfH,) 

l + # 

Ue dx (#+i) | 

[>-* 

CT 

i 

+ 

} 


dH 1 

3757 1 


1 

ffi(7 -1 )r A/2 R 2 ’ 

1 ' hr\ 

-H\ 1 

f 1 dm _ 1 C f\\ 
( 777 dx 2 0 ) J 


(7-1 )rM?R 2 ,, dH 

[ + /il 377Tj 



(9) 


( 10 ) 


9 


cICe 

dx 



R3 

H + H i 




1 + 0.075 Mi 


R i 


1 + 0.1 A/2 


( 11 ) 


where R\ = 1 T and R <2 = 1 + ^ De- 

fined in the following section, R$ is a factor that pro- 
vides transition between the equations for attached 
flow and separated flow. 

The subscript e in these equations refers to quan- 
tities at the boundary-layer edge, the subscript EQ 
denotes equilibrium conditions, and the subscript 
EQO denotes equilibrium conditions in the absence 
of secondary influences on the turbulent structure. 
(See ref. 24.) The parameters that appear in equa- 
tions (9) to (11) arc defined in the following section. 

Closure Conditions for Inverse Boundary 

Layer 

The inverse boundary-layer equations contain ad- 
ditional unknowns that must be specified by further 
assumptions (closure conditions) before the equa- 
tions can be solved. For separated flow, these clo- 
sure conditions are based upon the work of Melnik 
and Brook (ref. 17), which closely follows the analysis 


of LcBallcur in reference 15, where additional details 
may be found. 

The separated flow is represented by a detached 
free-shcar layer that is separated from the airfoil by a 
region of constant velocity reverse flow. The velocity 
profile (fig. 1) used to model the flow in the separated 
region is given by 

where 

c 2 = — ^ ( a i — 1; a 2 = i) 

ai (1 + 02*7*) V 2 ) 

6 * 

a = J 

Fp (r,) = 1 (0 <*?<*?*) 

F p (r,) = F r (fj) ( V * < r, < 1) 


4 


7) = 


7] = 


y-y 

1 -rj* 


The parameter a is determined by an iterative pro- 
cess described in a subsequent section of this paper. 
The function F c (rj) is Cole’s wake function: 


F C (v) = 2 (* + COS7r rf) 


The magnitude of the reversed flow is U m /U e — 
1 — C2. Following Melnik and Brook, 77* is given 

by 


rf (a) 


where 


f ba + (1 - b) 

(a m < a < 1) 

a (a - a s ) 2 

((ig < Oi < Om) 

1° 

(a < a s ) 

b = 2.1 


b 


e 

1 

r-H 

11 

e 

-5-1) 

a 5 = a\ 




The above assumptions, along with a mixing length 
formula for the Reynolds stresses, provide all infor- 
mation required for the following closure conditions 
(ref. 17): 


H (a) = {!-[! + B(rj*)] a}~ 1 


where 


B[rf) 


G 2 ff) 

G?0*) 


G 1 (»?*) = ai (1+02*7*) 
G 2 (/?*) = &i (1 + b 2 r 1 *) 



tfl(Z*) = tf(a)(l-l) 


( 12 ) 


(13) 


(C E ) EQO = 7T 2 0.0064C 2 

The derivative is calculated as follows: 

dH 

dH\ _ dH\ dH\ da _ dH\ OHi/da 

~ ~W + da W = ~W + Iwjdk 

The following equation for the skin friction Cj is the 
one used by Thomas in reference 25: 


(H < 2.4) 
(77 > 2.4) 


where 77* = 0 for attached flow. 

The following expression for the factor of 
equation (11) is used to provide a transition between 
the value used in reference 24 and the expression 
given in reference 17: 


r 1 74+0.31 H 


c f = 


(log i0 

+ (0.00011) [tanh (4 — 0^7^) - l] 


A = 


l 

1 - if 

1 1 

2 1 — if 


(On airfoil) 
(On wake) 


f 2.8 




0.15 


k (l-^) (0.8)(1 


(H<2) 
(H > 2) 


The expressions for the remaining parameters are 
identical to those given in the appendix. The equa- 
tions for F, C T} and (CV)eqo are slightly different 
from the corresponding equations in reference 17. 
The changes are based upon personal communica- 
tions with R. E. Melnik, the first author of refer- 
ence 17. As a result of these communications, the 
new expressions arc used in both the direct method 
and the inverse method. 


Viscous Boundary Conditions and Wake 

The coupling between the viscous boundary layer 
and the inviscid analysis is through the boundary 
conditions on the airfoil and wake. The boundary 
conditions given by equations (3) and (4) are modi- 
fied as follows: 


0? = S? + Sf+^ ± (x,e<x<x (e ;z = 0 ± ) (14) 

A<f> z = A (<5*) (z > x te ] z = 0*) (15) 

where 6* is the boundary-layer displacement thick- 
ness, the superscript ± refers to the upper or lower 
surface of the airfoil, and the A(. . .) denotes a jump 
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in the bracketed quantity across the wake. Equa- 
tion (14) does not include S* because of the quasi- 
steady assumption in the bound ary- layer equations. 

Numerical Implementation 

From the leading edge of the airfoil or wing, the 
boundary layer is approximated by the turbulent 
boundary layer on a flat plate. At a user-specified 
point, typically 10 percent chord, numerical inte- 
gration of the direct boundary-layer equations (6) 
to (8) is implemented with a fourth-order Runge- 
Kutta method. Downstream integration of the direct 
bound ary- layer equations continues until the flow 
nears separation, at which point the method switches 
to the inverse bound ary- layer equations (9) to (11). 
In the present application, the switch to the inverse 
equations occurs at H — 1.5. The inverse calcula- 
tion continues several chord lengths into the wake, 
even if the value of H drops below the switch value. 
The inverse equations can also be initiated at a user- 
specified point along the chord once integration of 
the direct equations has begun. 

At each time step, the inverse boundary-layer al- 
gorithm solves equation (12) by Newton iteration for 
a, given H. Then ff\ is determined from equa- 
tion (13) and the other parameters are computed. 
The CAP-TSD code also includes a subiteration ca- 
pability as part of the basic solution algorithm. With 
the boundary-layer calculations included, this sub- 
iteration results in successive viscous-inviscid itera- 
tions until the specified level of convergence has been 
achieved. 

Coupling between the inviscid outer flow solution 
and the direct boundary- layer calculation is straight- 
forward. Once the boundary-layer parameters are 
computed, the displacement thickness <5* required by 
the boundary conditions in equations (14) and (15) 
is given by 

<5* =6H 

For the inverse boundary-layer calculation, the dis- 
placement thickness <5* is computed by Carter’s 
method (ref. 20): 

^ew = Cl+^ld(^-l) 

where 

oj relaxation factor (typically 0.1 to 0.001) 

U ei inviscid velocity at boundary-layer edge 
Ue v viscous velocity at boundary-layer edge 


For the time-accurate calculations in this pa- 
per, values of u larger than 0.1 led to instabili- 
ties, although values larger than 1.0 are reported 
in reference 20, where only steady flow solutions are 
computed. 

Results and Discussion 

The present method has been applied to several 
test cases to evaluate its accuracy and range of ap- 
plicability. Some of these test cases, such as the 
NACA 64A010A and the NACA 0012 airfoil, have 
been calculated with previous codes (ref. 14) and 
are presented to confirm the accuracy of the present 
method as well as to demonstrate the improvements 
that have been obtained. The results for the buzz 
calculations for the airfoil on the P-80 aircraft repre- 
sent new applications of transonic small-disturbance 
theory with viscous-inviscid interaction. Figure 2 
contains profiles of the configurations studied. The 
NACA 64A010A airfoil has the coordinates of the 
section tested at the NASA Ames Research Center 
(ref. 26); this section had a small amount of camber 
and was slightly thicker than the symmetrical design 
section. 

Unless otherwise stated, the results for the two- 
dimensional calculations were obtained on a 142 x 84 
grid in x-z space. This grid extends ±20 chords in x 
and ±25 chords in z; it has 76 points on the airfoil. 
Also, the vorticity modeling option in the CAP-TSD 
computer code was turned on for all calculations 
except those for the airfoil on the P-80. 

NACA 64A010A Airfoil 

Ten AGARD computational test cases for the 
NACA 64A010A airfoil were calculated with a previ- 
ous version of the viscous-inviscid method 
(XTRAN2L) and compared with experimental re- 
sults in reference 14. In the present paper, the five 
cases that show the effect of frequency on unsteady 
lift and pitching-moment coefficients (i.e., cases 3 to 7 
listed in table I) are recalculated with the new com- 
puter code, and the results are compared with the 
previous calculations as well as with the experimen- 
tal results (ref. 26). The Mach number for these five 
cases was 0.796, the mean angle of attack was 0°, and 
the unsteady amplitude of harmonic oscillation was 
about 1°. The number of time steps per cycle used for 
the calculation was 720, the relaxation factor cj for 
the inverse calculations was 0.01, and the maximum 
number of subiterations was 20 with a convergence 
criterion of 0.0001. 

Figure 3 presents comparisons of the real and 
imaginary parts of the lift coefficient, and figure 4 
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presents similar comparisons of the pitching-moment 
coefficient (about the leading edge). In general, the 
CAP-TSD viscous results for the lift coefficient agree 
well with the experimental data. The real part of 
the lift coefficient is slightly overpredicted, and a 
small discrepancy exists in the imaginary part for 
intermediate frequencies. The imaginary part of the 
lift coefficient for the two lower frequencies (cases 3 
and 4) shows considerable improvement over the pre- 
vious calculations. An examination of the response 
time histories for the present results compared with 
those of previous calculations for cases 3 and 4 sug- 
gests that accurate calculation of the lift coefficient 
for these cases depends upon accurate calculation of 
a small separation bubble that develops at the base 
of the shock during the unsteady motion. 

As figure 4 shows, calculated results for the real 
part of the pitching-moment coefficient have the same 
trends with increasing frequency as the experimen- 
tal data, although the magnitudes are somewhat dif- 
ferent. The overprediction of the real part of the 
pitching-moment coefficient is consistent with the 
overprediction of the real part of the lift coefficient 
mentioned previously. The calculated imaginary part 
of the pitching-moment coefficient agrees fairly well 
with the experimental results with the largest differ- 
ences at the higher frequencies. 

NACA 0012 Airfoil 

The four AGARD cases for the NACA 0012 air- 
foil (table II) involve larger mean angles of attack (up 
to 4.86°) and larger amplitude pitch oscillations (up 
to 4.59°) than are normally considered appropriate 
for calculations with transonic small-disturbance the- 
ory. Calculated results for these cases are presented 
to investigate the range of applicability of the present 
theory. The grid for the NACA 0012 calculations is 
137 x 84 in x-z space, has 55 points on the airfoil, 
and extends ±20 chords in all directions. The relax- 
ation factor uj for the inverse calculations was 0.001, 
and the maximum number of subiterations was 20 
with a convergence criterion of 0.0001. The number 
of time steps per cycle used for the calculation was 
2048, and results were output every 32 time steps. 
For most time steps the convergence criterion was 
satisfied with 2 to 5 iterations, although all 20 itera- 
tions were usually required near the maximum angle 
of attack. 

Although viscous effects upstream of the shock 
are important physically, the interactive boundary 
layer tended to overpredict viscous effects when the 
boundary layer was initiated upstream of the shock. 
This overprediction can be the result of a laminar 


boundary layer upstream of the shock with tran- 
sition to turbulence at the shock. Because the 
present method assumes a turbulent boundary layer, 
the boundary layer was initiated downstream of the 
shock position. Because a change in shock position 
occurs during a cycle of motion, the boundary layer 
was initiated just downstream of the most rearward 
position of the shock during a cycle of oscillation. 
The exact location was determined by trial and error, 
and the inverse boundary layer was used from that 
point downstream. As subsequently shown, this ap- 
proach slightly underpredicts the viscous effects but 
yields overall results that agree well with the experi- 
mental data. 

Figures 5 to 9 compare the viscous calculations, 
inviscid calculations, and experimental data in refer- 
ence 27. Comparisons were made for instantaneous 
pressure distributions as well as for lift and moment 
coefficients versus angle of attack. Because of the 
difficulty in determining an appropriate time axis for 
the instantaneous pressure distribution comparisons 
during a harmonic oscillation, the experimental re- 
sults were Fourier analyzed to determine phase an- 
gles for each of the experimental points. The avail- 
able calculated results with the nearest phase angles 
were used for the comparisons. 

Figure 5 presents plots of the lift and pitching- 
moment coefficients versus a for the oscillatory 
cases 1, 2, 3, and 5. For cases 1 and 2, the vis- 
cous lift coefficient is slightly higher than the exper- 
imental data near the maximum angle of attack but 
agrees well with the experimental results elsewhere. 
For case 3, the viscous lift coefficient is slightly higher 
than the experimental results over the entire range 
of angles. For case 5, both the inviscid and viscous 
calculations give similar results. The lift coefficient 
is slightly lower than the experiment. Although the 
Mach number for case 5 is higher than those for the 
other cases, the mean angle of attack is lower and 
the shock is weaker. (See fig. 9.) Thus, the effect of 
the boundary layer is almost negligible. For cases 1 
to 3, the viscous lift and moment coefficients agree 
much better with the experimental results than do 
the inviscid calculations. 

The instantaneous pressure coefficients for 
cases 1, 2, 3, and 5 are compared in figures 6 to 9. 
In general, both the inviscid and viscous calculations 
compare well with the experimental data. The most 
noticeable difference between the calculated values 
and the experimental results is the overprediction 
of the leading-edge suction pressure by both invis- 
cid and viscous calculations. This overprediction of 
leading-edge pressure was not present in previous cal- 
culations with the inviscid code by Batina (ref. 19). 
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The source of the discrepancy is not known, although 
it may result from the different grids used in the 
calculations. During the quarter cycle after maxi- 
mum a (figs. 6(c), 7(c), and 8(d)), the shock posi- 
tion for the viscous calculation is noticeably better 
than that for the inviscid calculation as a result of 
separation effects. Otherwise, the viscous and invis- 
cid results arc comparable. For case 5, the viscous 
and inviscid results are essentially the same. 

The results presented herein for the NACA 0012 
airfoil demonstrate that transonic small-disturbance 
theory with an interactive inverse boundary layer can 
predict with reasonable accuracy the air loads due to 
moderately large-amplitude pitch oscillations. 

P-80 Aileron Buzz 

Flight test measurements on the P-80 fighter that 
were conducted during the mid 1940’s demonstrated 
that a limit cycle oscillation of the aileron control 
surfaces occurred during transonic flight conditions 
(ref. 28). This limit cycle oscillation of control sur- 
faces is commonly referred to as aileron buzz. Wind 
tunnel measurements of a partial span P-80 wing 
were performed in the NASA Ames 16-Foot High- 
Speed Wind funnel and demonstrated that the ba- 
sic physical mechanism driving the oscillation is a 
lag in hinge moment that follows control surface dis- 
placement (ref. 29). A two-dimensional calculation 
that used the P-80 airfoil section shown in figure 2, 
an NACA 65 i -2 13 ( a = 0.5) airfoil, was performed in 
reference 3 with a Navier-Stokes method for the aero- 
dynamic loads. Although the calculations by Steger 
and Bailey were exploratory, they demonstrated that 
aileron buzz can be studied with Navier-Stokes cal- 
culations. In the present paper, numerical results 
demonstrate that aileron buzz can also be investi- 
gated through use of transonic small-disturbance the- 
ory with an interactive inverse boundary layer. The 
aileron is modeled as a single control surface mode 
shape pitching about the three-quarter-chord loca- 
tion. The physical quantities that define the model 
are taken from reference 28: Moment of inertia = 
0.4083 ft-lb/sec 2 , Mean chord = 4.83 ft, and Aileron 
span = 7.5 ft. 

The P-80 calculations that were performed in- 
cluded the effects of shock generated entropy, and the 
inverse boundary- layer calculation was initiated at 
12 percent chord. A weak spring was inserted at the 
aileron hinge because the computer code does not al- 
low zero spring stiffness. The Reynolds number used 
in the calculations was 20.0 million. No subiterations 
were used during the calculations and the relaxation 
factor u; was 0.1. Results were calculated for three 
values of oq(— 1°, 0°, +1°), and the time step was 


varied from 214 to 825 steps per cycle as subsequently 
discussed. To determine aileron buzz, a steady so- 
lution was first calculated at a Mach number close 
to the anticipated buzz condition. This steady solu- 
tion was then used as a starting solution for an aero- 
elastic calculation with the dynamic pressure q fixed. 
The Mach number was varied until a buzz oscillation 
was obtained with the aileron released from its unde- 
flccted position. This procedure does not necessarily 
determine the minimum Mach number for the onset 
of aileron buzz. According to reference 3, buzz can 
be induced by releasing the aileron from a deflected 
position at conditions where it did not buzz when 
released from the undeflectcd position. 

The buzz calculations were found to be highly 
nonlinear. Small changes in the parameters of the 
problem can significantly affect the numerical results. 
The nonlinear variation of static hinge moment with 
Mach number is undoubtedly responsible for much 
of the nonlinearity. As noted in reference 29, in the 
transonic speed range the “hinge moments are a sen- 
sitive function of Mach number.” The present anal- 
ysis confirms this observation. Figure 10 shows the 
calculated aileron hinge moment versus Mach num- 
ber for three airfoil angles of attack. These results 
are for steady flow conditions with the aileron held 
fixed at its undeflected position. As figure 10 shows, 
for Mach numbers higher than 0.78, the aileron hinge 
moment is highly nonlinear. This strong nonlinearity 
results in unsteady flow solutions that are sensitive 
to some of the parameters used to obtain the numer- 
ical results. This effect is discussed further in the 
following paragraphs. 

Figure 11 shows the steady pressure distributions 
of the upper and lower surfaces for M = 0.82 and 
ol = 0°. The shock on the upper surface is at 78 per- 
cent chord, just downstream of the aileron hinge line. 
The lower surface shock is at 67 percent chord. The 
difference between the upper and lower surface pres- 
sures over the aileron results in a small moment about 
the aileron hinge. When the aileron is released at the 
start of an aeroelastic calculation, this unbalanced 
hinge moment deflects the aileron upwards and an 
unsteady oscillation begins. When the Mach number 
is increased slightly to M = 0.8243, this unsteady os- 
cillation develops into a buzz limit cycle. Figure 12 
shows the unsteady aileron deflection angle /3 versus 
time for M = 0.8243 during an aeroelastic calcula- 
tion. As the figure shows, after three cycles of oscil- 
lation the aileron response settles into a limit cycle 
oscillation. The reduced frequency is k = 0.3808, 
which corresponds to a frequency of 21.67 Hz. This 
value compares well with the wind tunnel values 
that varied between 19.4 Hz and 21.2 Hz over a 
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variety of test conditions (ref. 29). The buzz fre- 
quency reported during the flight tests was 28 Hz 
(ref. 28). The calculated amplitude for the limit cycle 
oscillation is about ±2° about an aileron uplift angle 
of -1°. In the wind tunnel tests and the Navier- 
Stokes calculations (refs. 29 and 3), the unsteady 
amplitude is about ±10°, and the value reported in 
the flight tests is 2° (ref. 28). This buzz calculation 
appears to represent the current limit of the inverse 
boundary-layer code as the calculation diverges when 
the Mach number is slightly increased. 

Also, the calculated buzz conditions for this ex- 
ample changed slightly with the value of the time 
step used in the numerical integration. A summary 
of the variation observed for this example is pre- 
sented in table III. The buzz Mach number varied 
from M = 0.8241 for At = 0.02 to M = 0.8247 
for At = 0.04. The buzz frequency varied between 
20.9 Hz for At = 0.04 and 21.7 Hz for At = 0.01. Al- 
though fully converged solutions were not obtained 
for this highly nonlinear case, the buzz phenomena 
persisted as the time step was decreased. 

Figure 13 presents comparisons of buzz bound- 
aries calculated by this method with the Navier- 
Stokes calculations of Steger and Bailey (ref. 3) and 
the wind tunnel experiments of reference 30. The 
CAP-TSD results were obtained for two values of 
the dynamic pressure: q « 644 psf, corresponding 
to sea level at standard atmospheric conditions, and 
q « 293 psf, corresponding to an altitude of 6096 m 
(20000 ft). The time step for the CAP-TSD cal- 
culations was = 0.02. As mentioned previously, 
the CAP-TSD results were obtained by releasing the 
aileron from its undeflected position and do not nec- 
essarily represent the minimum Mach number for the 
onset of buzz. For q « 644 psf, the buzz frequency 
was about 22 Hz for all cases. For q » 293 psf, 
the calculated buzz frequency was reduced to about 
13 Hz. The Navier-Stokes result at a = —1° and 
M — 0.83 was obtained by releasing the aileron from 
its undeflccted position, and this result is in reason- 
able agreement with the present calculations. 

A CAP-TSD calculation that included three cy- 
cles of acroclastic oscillations required less than 
10 minutes of computer time on the Cray-2 com- 
puter at NASA Langley Research Center. Thus, 
the present technique offers an economical and rea- 
sonably accurate method for studying aileron buzz 
oscillations. 

F-5 Wing 

The CAP-TSD code with a stripwise viscous 
boundary layer has been applied to both steady and 


unsteady cases for the F-5 wing. The F-5 wing has a 
panel aspect ratio of 1.58, a leading-edge sweep angle 
of 31.9°, and a taper ratio of 0.28. The airfoil for this 
wing is a modified NACA 65A004.8 airfoil that has 
a drooped nose and is symmetrical aft of 40 percent 
chord. (See fig. 2(b).) The grid used for the F-5 cal- 
culations contains 137 x 30 x 84 points in the x, y, 
and z directions. This grid extends ±20 chords in the 
x and z directions and 2 semispans in the y direction, 
and it has 20 stations along the wing semispan with 
55 points on each chord. The experimental data used 
in the comparisons arc from reference 30. 

Figure 14 presents comparisons of experimental 
and calculated pressure distributions for steady flow 
at a Mach number of 0.897 and c\ = ^0.004°. The 
inviscid calculations for the two inboard stations 
(q s = 0.181 and 0.355) indicate a mild shock on the 
upper surface, whereas both the viscous calculation 
and the experiment show no evidence of a shock at 
these stations. The viscous calculation indicates the 
development of a shock at q s = 0.512, and the ex- 
perimental data suggest a mild shock at tj s = 0.641. 
All three results show a shock at the four outboard 
spanwise stations. The viscous shock is located about 
2 percent upstream of the inviscid shock and gen- 
erally is in better agreement with the experimen- 
tal data, although the lack of experimental points 
makes the experimental location of the shock uncer- 
tain. Near the wing tip ( r] s — 0.977), both the viscous 
and the inviscid calculations predict a shock location 
slightly upstream of the experimental results. The 
differences between calculated and experimental re- 
sults near the wing tip may result from slight differ- 
ences between the analytical model and the experi- 
mental wing in this region, highly three-dimensional 
flow effects near the wing tip, or coarseness of the 
grid used for the calculations. 

Unsteady calculations were made for a Mach 
number of 0.899, qq = 0.002°, a i = 0.109°, and 
k = 0.137 (/ — 20 Hz). (See figs. 15 and 16.) 
The time step used in the calculations corresponds 
to 500 steps per cycle of motion, and five cycles were 
calculated to allow for decay of any initial transients. 
The last cycle of the calculated results was Fourier 
analyzed to determine the harmonic content, and the 
first-harmonic components are compared with the 
experimental data in figures 15 and 16. The most 
noticeable feature of the upper surface unsteady 
pressures in figure 15 is the large variation in the 
calculated shock pulses near midchord. Although the 
viscous results are not always closer to the experi- 
mental data points in this region, the maximum am- 
plitudes of the viscous results arc substantially less 
than those of the inviscid results, and the viscous 
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shock pulses are slightly upstream of their inviscid 
counterparts. Hence, the viscous calculations are in 
better agreement with the data than the inviscid re- 
sults. As also shown in figure 15, viscosity has little 
effect on the upper surface pressures upstream of the 
shock wave. Immediately downstream of the shock 
wave, the viscous results generally agree better with 
the experimental data than do the inviscid results. 
Near the trailing edge, essentially no difference occurs 
between the calculated viscous and inviscid unsteady 
pressures on the upper surface of the wing. Figure 16 
indicates that on the lower surface of the wing the 
leading-edge suction peak is poorly predicted by both 
inviscid and viscous calculations. Away from the 
leading edge, the only significant differences between 
the viscous and inviscid unsteady calculations occur 
in the vicinity of midchord. In this region, the viscous 
results agree better with the experimental data than 
do the inviscid calculations. For the F-5 wing, the 
viscous results calculated with the interactive strip- 
wise boundary layer provide a qualitative indication 
of the effects of viscosity within the cost-effective 
framework of transonic small-disturbance theory. 

Conclusions 

A method is presented for calculating turbu- 
lent viscous effects for two- and three-dimensional 
unsteady transonic flows, including flows involving 
mild separation. The method uses Green’s lag- 
entrainment equations for attached turbulent flows 
and an inverse boundary- layer method developed by 
Melnik and Brook for separated turbulent flows. The 
inverse boundary-layer equations are coupled with 
the inviscid flow calculation through use of Carter’s 
method. The viscous method uses steady boundary- 
layer equations in a quasi-steady manner, and the 
three-dimensional viscous effects are included in a 
stripwise fashion. The method has been applied to 
several two-dimensional test cases as well as a three- 
dimensional wing planform. The applications include 
the calculation of limit cycle oscillations, such as 
those that occur in aileron buzz. Comparisons are 
presented with experimental data and inviscid anal- 
yses as well as with another interactive boundary- 
layer method and Navier- Stok es calculations. The 
results demonstrate that accurate solutions are ob- 
tained for unsteady two-dimensional transonic flows 
with mild separation, and qualitative viscous effects 
are predicted for three-dimenisonal flow fields. The 
results have led to the following general conclusions: 

T For the NACA 64A010A airfoil, the lift coefficient 
calculated with the CAP-TSD viscous code shows 


a considerable improvement over previous calcu- 
lations for the lower frequency AGARD cases. 
In particular, for the two lower frequencies, the 
imaginary part of the lift coefficient calculated 
with the CAP-TSD viscous code agrees well with 
the experimental data. 

2. For the NACA 0012 airfoil, reasonably accurate 
calculations of lift and moment coefficients have 
been obtained with the viscous code for large- 
amplitude pitch oscillations, which are usually 
considered outside the range of transonic small- 
disturbance theory. 

3. Instantaneous shock positions calculated with the 
viscous code during large-amplitude pitch oscilla- 
tions of the NACA 0012 airfoil show better agree- 
ment with the experimental results than do those 
of the inviscid code. 

4. The CAP-TSD viscous code accurately calcu- 
lated the Mach number and frequency for the on- 
set of aileron buzz for the airfoil on the P-80. 
The buzz boundary calculated with the viscous 
code agrees well with the experimental data and 
Navier-Stokes calculations. 

5. The viscous code is relatively economical in terms 
of computer time. For example, an aeroelastic 
buzz calculation for three cycles of motion re- 
quired less than 10 minutes computer time on a 
Cray-2. 

6. For the F-5 wing, steady calculations with the vis- 
cous code predicted shock locations and strengths 
in better agreement with experimental results 
than did the inviscid calculations. Unsteady 
calculations indicate that the stripwise viscous 
boundary layer can provide a qualitative indi- 
cation of viscous effects within the cost-effective 
framework of transonic small-disturbance theory. 

7. The results presented demonstrate that the vis- 
cous solutions computed with the present algo- 
rithm can provide predictions of pressure dis- 
tributions for unsteady transonic flow involving 
moderate-strength shock waves and mild flow sep- 
aration that correlate better, sometimes signifi- 
cantly better, with experimental values than do 
the inviscid solutions. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
April 20, 1992 
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Appendix 

Viscous Parameters for Attached Flow 


For attached boundary layers, the various depen- 
dent variables and functions are evaluated from the 
following expressions. 
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The additional parameters required to specify the 
boundary-layer equations completely, together with 
the default values in the code (in parentheses), arc 

the free-stream chord Reynolds number iVp e ^10 7 ^, 

the free-stream temperature T in kelvins (300 K), 
the turbulent Prandtl number JVp r ^ (0.9 for air), and 
the Sutherland law viscosity constant N$ u in kelvins 
(110 K for air). 
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Table I. Analytical Test Cases for NACA 64A010A Airfoil 
[M = 0.796; N Re = 12.5 x 10 6 ; a 0 = 0°; x a /c = 0.25] 


<*1, deg 

/, Hz 

k 

1.03 

4.2 

0.025 

1.02 

8.6 

.051 

1.02 

17.2 

.101 

1.01 

34.4 

.202 

.99 

51.5 

.303 


Table II. Analytical Test Cases for NACA 0012 Airfoil 
[k = 0.081; x„/c = 0.25] 


U, m/sec 


A r R, 

4.8 x 10 6 
4.8 
4.8 
5.5 


f*0i deg 


£*l, deg 


Table III. Variation of Buzz Conditions With Time Step at a = 


Steps per cycle Mg Qn- psi 

2l4 0.8247 644.9 

422 .8241 644.3 

825 -8243 644.5 
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Figure 6. Comparison of instantaneous unsteady pressures for case 1 for the NACA 0012 airfoil. M = 0.60 
qo = 2.89°; r*i = 2.41°; k = 0.081; A r Re - 4.8 x 10 6 . 
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Figure 7. Comparison of instantaneous unsteady pressures for case 2 for the NACA 0012 airfod. M - 0.59 
qo = 3.16°; «i = 4.59°; k = 0.081; N Re = 4.8 x 10 6 . 
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Figure 8. Comparison of instantaneous unsteady pressures for case 3 for the NACA 0012 airfoil. M — 0.5S 
rvn = 4.86°; «i = 2.44°; k = 0.081; A r Re = 4.8 x 10 6 . 
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Figure 9. Comparison of instantaneous unsteady pressures for case 5 for the NACA 0012 airfoil. M — 0.75 
an = 0.02°; <*i = 2.51°; k = 0.081; N Ke = 5.5 x 10 6 . 
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Figure 10, Variation of calculated aileron hinge moment with Mach number for airfoil on the P-80 at three 
angles of attack with undeflected flap. 



Figure 11. Steady pressure distribution for airfoil on the P-80. M — 0.82; cto — 0°. 



Figure 12. Typical calculated aileron buzz oscillation for airfoil on the P-80 released from undeflected position. 
M = 0.8243; a 0 = 0°; At = 0.01. 
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Figure 14. Comparison of calculated and experimental : 
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; pressure distributions for F-5 wing. M = 0.89 
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Figure 15. Comparison of calculated and experimental unsteady pressure distributions for upper surface 
of F-5 wing. M = 0.899; k = 0.137; q 0 - 0.002°; <*1 = 0.109°. 
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Figure 15. 
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Figure 16. Comparison of calculated and experimental unsteady pressure distributions for lower surface of F-5 
wing. M = 0.899; k = 0.137; ag = 0.002°; ai = 0.109°. 
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Figure 16. Concluded, 
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